A deeper understanding of system interactions can explain contradictory field results on pesticide impact on honey bees

While there is widespread concern regarding the impact of pesticides on honey bees, well-replicated field experiments, to date, have failed to provide clear insights on pesticide effects. Here, we adopt a systems biology approach to gain insights into the web of interactions amongst the factors influencing honey bee health. We put the focus on the properties of the system that depend upon its architecture and not on the strength, often unknown, of each single interaction. Then we test in vivo, on caged honey bees, the predictions derived from this modelling analysis. We show that the impact of toxic compounds on honey bee health can be shaped by the concurrent stressors affecting bees. We demonstrate that the immune-suppressive capacity of the widespread pathogen of bees, deformed wing virus, can introduce a critical positive feed-back loop in the system causing bistability, i.e., two stable equilibria. Therefore, honey bees under similar initial conditions can experience different consequences when exposed to the same stressor, including prolonged survival or premature death. The latter can generate an increased vulnerability of the hive to dwindling and collapse. Our conclusions reconcile contrasting field-testing outcomes and have important implications for the application of field studies to complex systems.

In particular, according to Theorem 5, bistability in terms of high and low bee health can arise in a given interval of immune suppression capacity, below which the system is monostable at high bee health and above which the system is monostable at low bee health. In the range of bistability, the final state depends on the system initial conditions and minimal differences in these conditions can decide whether the final state is at high or low bee health.

Supplementary Methods
First to proceed, we note that Fig. 2 As it can be deduced from the proof of Proposition 1 in the section "Structural analysis of the bee health model" of the main text, system (1) is monotone (since the matrix M is Metzler 1 ).
Monotone systems can exhibit multiple equilibria and hence multistability 1 . To investigate whether this is the case for system (1) we must specify the form of the functions g, f , f and h. The following choices are among the simplest capable of portraying the desired features (monotonicity, convexity, saturation, etc.). As such, they reasonably represent the concerned interactions without compromising the generality of the conceptual model. We thus focus on where we set δ = 1 for each state variable without loss of generality.
The search for equilibria passes through annihilating the right-hand side, leading to the algebraic system A full analysis of (3) is not immediate. We thus proceed by analyzing different subsystems, finally showing that the presence of pathogens impairing the immune system is a necessary condition for multistability. In the sequel, with the term equilibrium we implicitly refer to positive values of the state variables considered in each specific subsystem, the only biologically meaningful. In this sense, all the following results share the common assumption of having u T sufficiently small. This assumption prevents the honey bee health at equilibrium from being null or negative (and the same for the parasites when present). We stress that the assumption is necessary because a large u T would be disruptive and it would prevent any survival possibility.
We start by considering the case in which the only stress factor is a toxic compound.
Proposition 2: If u T is sufficiently small, the subsystem of honey bees and toxic compounds (i.e., first and second equation in (2) with ν = ρ = c 2 = 0) has a unique equilibrium, globally asymptotically stable.
It is immediate to recover the nullcline † for x HB and the nullcline for x T C Equating the two leads to finding that possible equilibria are characterized by x T C = y for y any positive solution of where LHS(y) := ξ − ωy 1 ω + µ + σ − c 6 u T (1 + βy)(1 + c 3 y)(1 + c 5 y) † A nullcline for a state variable of a planar system is a curve in the phase plane along which that variable is locally constant. As such it is obtained by annihilating the relevant derivative and hence the right-hand side of the corresponding ODE. It follows that equilibria are intersection points of nullclines of different state variables. and RHS(y) := ωy α(1 + c 3 y)(1 + c 5 y) + c 0 u S (1 + c 1 u C )(1 + βy)(1 + c 5 y) On the one hand, LHS is a fourth degree polynomial in y, with negative leading term, three negative roots and one positive root .
The fact that y + is positive follows by assuming u T < (µ + σ)/c 6 , condition guaranteeing that the equilibrium coordinate x HB is positive (indeed, the minimum such value is µ + σ − c 6 u T when x T C → +∞ as it follows from (5)). On the other hand, RHS is a third degree polynomial in y, with positive leading term and a root in 0. By factoring this trivial root out, the remaining second degree polynomial has positive coefficients, hence no positive roots (as it follows by applying Descartes' rule of signs ‡ ). It follows that there exists a unique positive solution y of (7), being indeed LHS(0) = ξ > 0 = RHS(0).
Global stability of the corresponding equilibrium immediately follows by analyzing the nullclines. In the phase plane (x T C , x HB ), orbits through points above (5) move downward (sinceẋ HB < 0), while orbits through points below move upward (sinceẋ HB > 0). On the other hand, orbits through points to the right of (6) move leftward (sinceẋ T C < 0), while orbits through points to the left move rightward (sinceẋ T C > 0). Their unique intersection point (i.e., the equilibrium) is thus globally attracting, Supplementary Fig. 2a. This ends the proof of Proposition 2.
Going back to the general model (1), from the proof above we can conclude that any reasonable choice of the functions g, f and f satisfying the required assumptions of smoothness, monotonicity and convexity can determine uniqueness and global stability. ‡ The number of positive roots of a polynomial is at most the number of sign changes in the series of ordered coefficients in its canonical expression.
The second result concerns the case in which the only stress factor is due to the presence of parasites.
Proof: Similarly to the proof of Proposition 2 we reduce (2) to and recover the nullcline for x HB and the nullcline for Equating the two leads to finding that possible equilibria are characterized by x V A = z for z any positive solution of where and On the one hand, LHS is a third degree polynomial in z, with negative leading term, two negative roots and one positive root In fact, from (9) we have x HB ∈ (x HB,min , x HB,max ) for Then, by assuming u T sufficiently small (in particular, smaller with the last inequality holding true, otherwise x V A would not be positive from (10). On the other hand, RHS is a second degree polynomial in z, with positive leading term and two negative roots. Since holds true for u T sufficiently small, it follows that there exists a unique positive solution z of (11).
Finally, that the corresponding equilibrium is globally asymptotically stable follows similarly as in the proof of Proposition 2, Supplementary Fig. 2b. This ends the proof of Proposition

3.
A comment similar to the one placed soon below the proof of Proposition 2 holds here as well: any reasonable choice of the functions g, h, f and f in the first and third equation of the general model (1) satisfying the required assumptions can ensure uniqueness and global stability.
The third result concerns the case in which the stress factor is due to the presence of both toxic compounds and parasites. As the relevant subsystem is not planar anymore, the proof requires slightly different and more geometrical considerations.
Proposition 4: If u T is sufficiently small, the subsystem of honey bees, toxic compounds and parasites (i.e., first, second and third equation in (2) with ρ = θ = 0) has a unique equilibrium, globally asymptotically stable.
and note that u T must be sufficiently small similarly to the proofs of Proposition 2 and Proposition 3. Recover x HB from the last equation as In order to prove uniqueness, in the phase space (x T C , x V A , x HB ) we first consider from the first of (12) the family of curves curve is smooth, bounded, positive and convex, monotonically decreasing from and both these quantities are monotonically decreasing w.r.t. x V A , Supplementary Fig. 3a. We consider as well the family of curves (13): each curve is a line of slope and is monotonically increasing from and both these quantities are again monotonically decreasing w.r.t. x V A , Supplementary Fig. 3b.
It follows that for any fixed positive value of x V A , there is a unique intersection point between Now, we know from Proposition 3 that in the plane (x V A , x HB ) (i.e., x T C = 0) the curves have a unique intersection with positive coordinate x V A . Beyond this intersection g 1 is above g 2 and the former decreases more slowly than the latter, so that the above mentioned unique intersection point between f 1 and f 2 decreases with x V A , Supplementary Fig. 3c.
Finally, by following arguments similar to those used in the proof of Proposition 2 it is not difficult to show that, when projected into the plane (x T C , x HB ), the locus of such points forms a decreasing curve which intersects the curve given by the second of (12) only once, thus proving uniqueness, Supplementary Fig. 3d. where is introduced for brevity.
Theorem 5: Let u T be sufficiently small and a(ε) and b(ε) be given by (14). If then, by varying ε ≥ 0, the subsystem of honey bees and pathogens (i.e., first and fourth equation in (2) with β = ν = χ = c 2 = c 3 = c 5 = 0) exhibits two fold bifurcations § , viz. at ε = ε min and ε = ε max with 0 < ε min < ε max , and (A1) a unique equilibrium, globally asymptotically stable with high bee health, if and only if ε < ε min ; (A2) three equilibria, two locally asymptotically stable (with high and low bee health, respectively) and one unstable (with intermediate bee health), if and only if ε min < ε < ε max ; (A3) a unique equilibrium, globally asymptotically stable with low bee health, if and only if ε > ε max . § A fold bifurcation occurs when there exists a value of a system parameter below which, locally, there are no equilibria and above which, always locally, there are two equilibria (which coincide at the bifurcation value), or the other way around 2 .
Otherwise, if (15) does not hold, the system always has a unique equilibrium, globally asymptotically stable.
Proof: Again, the necessity of the assumption that u T is sufficiently small is clear from the proofs of Proposition 2, Proposition 3 and Proposition 4. In the usual manner we recover the nullcline for x HB and the nullcline for Equating the two leads to finding that possible equilibria are characterized by x V I = v for v any positive root of the third-degree polynomial parametrized by ε.
We start from the case ε = 0, by observing that p(v; 0) is a parabola heading upward and p(0; 0) = −γ < 0 implies the existence of one positive root. This root is unique since a(0) > 0, so that p(v; 0) has only one sign change independently of the sign of b(0) and the conclusion follows by applying Descartes' rule of signs.
Assume then ε > 0. We show that b(ε) ≤ 0 implies uniqueness (the same reasoning is valid also for a(ε) ≥ 0, the proof of which we thus omit). Indeed, in this case p(v; ε) has only one sign change (independently of the sign of a(ε)), hence at most one positive root again according to Descartes' rule of signs, and it has exactly one positive root since p(0; ε) = −γ ≤ 0 and the leading term is positive.
So b(ε) > 0 is a necessary condition to violate uniqueness, which we assume to hold true in the following. By considering only simple roots for the time being (i.e., those for which p (v; ε) = 0), uniqueness is lost only in favor of three positive roots. For this to be the case there must be three sign changes in the coefficients of p(v; ε), which is possible only if a(ε) < 0.
Therefore this constraint is a second necessary condition for non-uniqueness. A third necessary condition is given by the existence of a local minimum and a local maximum, which is equivalent to p (v; ε) = 3ρεv 2 + 2a(ε)v + b(ε) having two distinct zeros. Hence it must be b(ε) < a 2 (ε)/(3ρε), thus completing (15).
The final result under the last condition follows by observing first that p(γ; ε), p (γ; ε) and p (γ; ε) are all positive, so that the rightmost positive root of p(v; ε) as well as its local minimum and maximum lie all to the left of γ and second that when ργ > 1 (which is necessarily implied by a(ε) < 0) p(v; ε) decreases with ε when v < γ, in fact this implies The outcome is summarized in Fig. 2d, which illustrates the behavior of p(v; ε) under (15) for increasing ε (i.e., from top to bottom ε = 1.55: one equilibrium at high bee health; ε = ε min ≈ 1.56519: one equilibrium at high bee health and fold bifurcation at low bee health; ε = 1.58: three equilibria; ε = ε max ≈ 1.61591: one equilibrium at low bee health and fold bifurcation at high bee health; ε = 1.63: one equilibrium at low bee health; see Supplementary Tab. 4 for the values of the parameters), and fully justifies the thesis with the understanding that a low level of x V I (i.e., a low value of v) corresponds to a high level of bee health. In particular, when ε = 0 the only positive root is to the left of the local maximum (low v, high bee health).
When ε increases a first value is encountered, i.e., ε min , at which the local minimum touches the v-axis from above, generating two new positive roots to the right of the existing one (hence with intermediate and low bee health). By increasing ε further these three positive roots continue to exist, with the intermediate one moving leftward to the leftmost one, until at a second value for ε, i.e., ε max , these two roots collide and coincide with the local maximum to then disappear, leaving alone the positive rightmost root (high v, low health). This ends the proof of Theorem

5.
As noted in the proof above, a(ε) < 0 cannot hold true if ργ ≤ 1, which is biologically meaningful: when the effect of the pathogen's growth γ combined with the influence ρ on the bee health is sufficiently low, the immune suppression capacity is not effective. Moreover, when ργ is large enough so that a(ε) < 0 it is not difficult to show that there are choices of the parameter values for which also the second condition in (15) is satisfied, so that multiple equilibria are indeed possible (as in Fig. 2c). In particular, for this to be true also σ must be large enough, which in fact weighs the negative direct effect of pathogens on bee health (through ρ).
To complete the whole analysis, we observe that also the full system (1) with ε = 0 exhibits uniqueness and global stability, as the interaction of pathogens and parasites with the bee are modeled in the same way in this case. Instead, if ε > 0, multistability in the sense of Theorem 5 is possible, as it can be proved by assessing the introduction of toxic compounds and/or parasites on the base of continuity arguments. The mite is adapted to the optimal temperature for honey bees and may suffer from non-optimal temperatures; however, it is extremely rare that brood -where the mite reproduces-can suffer sub-optimal temperatures 25 Supplementary Table 3: Properties of the functions of the model of honey bee health, together with a summary of the biological effects they account for and a reference to the conceptual model in Fig. 1.  Fig. 1 g smooth, bounded, positive, convex and decreasing to 0    Supplementary Fig. 2b Supplementary Fig. 3 τ HB 1 1